/*
** ComplexFloatSymmetricLevinson.cs
**
** Class that implements UDL factorisation of the inverse of a symmetric square
** Toeplitz matrix using the Levinson Algorithm.
**
** Copyright (c) 2003-2004, dnAnalytics Project. All rights reserved.
*/

#region Using directives

using System;
using dnA.Exceptions;

#endregion Using directives

namespace dnA.Math
{

    /// <summary>
    /// A Levinson solver for a symmetric square Toeplitz systems.
    /// </summary>
    /// <threadsafety static="true" instance="true"/>
    /// <remarks>
    /// This class implements an <B>UDL</B> decomposition of the inverse of a symmetric 
    /// square Toeplitz matrix. The decomposition is based upon Levinson's algorithm. As
    /// a consequence all operations require approximately <B>N</B>*<B>N</B> FLOPS, where
    /// <B>N</B> is the matrix order. This is significantly faster than Cholesky's
    /// factorization for symmetric matrices.
    /// <para>
    /// One disadvantage of Levinson's algorithm is that all the principal sub-matrices
    /// of the Toeplitz matrix must be non-singular. During the decomposition,
    /// sub-matrices are checked to ensure that they are non-singular. If one of the
    /// sub-matrices is near-singular, then the accuracy of the decomposition will
    /// be degraded.
    /// </para>
    /// <para>
    /// A outline of this approach to the UDL decomposition of inverse Toeplitz matrices is
    /// found in the following reference:
    /// </para>
    /// <para>
    /// <EM>Sun-Yuan Kung and Yu Hen Hu</EM>, A Highly Concurrent Algorithm and Pipelined
    /// Architecture for Solving Toeplitz Systems, IEEE Transactions on Acoustics,
    /// Speech and Signal Processing, Volume ASSP-31, Number 1, Febuary 1983, pages
    /// 66 - 75.
    /// </para>
    /// </remarks>
	sealed public class ComplexFloatSymmetricLevinson : Algorithm
    {

        #region Fields

		/// <summary>
		/// The left-most column of the Toeplitz matrix.
		/// </summary>
		private readonly ComplexFloatVector m_ToeplitzColumn;

		/// <summary>
		/// The system order.
		/// </summary>
		private readonly int m_Order;

		/// <summary>
		/// Flag to indicate if system is singular.
		/// </summary>
		private bool m_IsSingular;

		/// <summary>
		/// lower triangular matrix.
		/// </summary>
		private ComplexFloat[][] m_LowerTriangle;

		/// <summary>
		/// The diagonal vector.
		/// </summary>
		private ComplexFloat[] m_Diagonal;

        #endregion Fields

        #region Properties
		
		/// <summary>
		/// Check if the Toeplitz matrix is singular.
		/// </summary>
		/// <remarks>
		/// This property actually tells us if any of the principal sub-matrices are
		/// singular as well as the main matrix.
		/// </remarks>
        public bool IsSingular
        {
            get
            {
				Compute();
				return m_IsSingular;
            }
        }

        #endregion Properties

        #region Constructors

		/// <summary>
		/// Qualified constructor with specified Toeplitz vector.
		/// </summary>
		/// <param name="T">The left-most column of the Toeplitz matrix.</param>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> is zero.
		/// </exception>
        public ComplexFloatSymmetricLevinson(ComplexFloatVector T)
        {
            // check parameter
            if (T == null)
            {
                throw new System.ArgumentNullException("T");
            }
            else if (T.Length == 0)
            {
				throw new RankException("The length of T is zero.");
            }

			// save the vector
			m_ToeplitzColumn = T.Clone();
			m_Order = m_ToeplitzColumn.Length;

			// allocate memory for lower triangular matrix
			m_LowerTriangle = new ComplexFloat[m_Order][];
			for (int i = 0; i < m_Order; i++)
			{
				m_LowerTriangle[i] = new ComplexFloat[i+1];
			}

			// allocate memory for diagonal
			m_Diagonal = new ComplexFloat[m_Order];

        }


        #endregion Constructors

		#region Protected Members

		/// <summary>
		/// Calculate the UDL decomposition of the inverse Toeplitz matrix.
		/// </summary>
		/// <remarks>
		/// The member checks each iteration for a singular principal sub-matrix.
		/// If a singular sub-matrix is detected, the private boolean <EM>m_IsSingular</EM>
		/// is set to <B>true</B> and the decomposition halted.
		/// <para>
		/// When the decomposition is completed, the lower triangular matrix is stored
		/// in the jagged array <EM>m_LowerTriangle</EM>, the prediction coefficients
		/// in the array <EM>m_Diagonal</EM> and the private boolean <EM>m_IsSingular</EM>
		/// is set to <B>false</B>.
		/// </para>
		/// </remarks>
		protected override void InternalCompute()
		{

			int i, j, l;			// index/loop variables
			ComplexFloat Inner;			// inner product
			ComplexFloat K;				// reflection coefficient
			ComplexFloat[] B;				// reference to previous order coefficients
			ComplexFloat[] A;				// reference to current order coefficients

			// check if principal diagonal is zero
			if (m_ToeplitzColumn[0] == ComplexFloat.Zero)
			{
				m_IsSingular = true;
				return;
			}

			// setup zero order solution
			B = m_LowerTriangle[0];
			B[0] = ComplexFloat.One;
			m_Diagonal[0] = m_ToeplitzColumn[0];

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++, B = A)
			{
				// calculate inner product
				Inner = ComplexFloat.Zero;
				for ( j = 0, l = 1; j < i; j++, l++ )
				{
					Inner += m_ToeplitzColumn[l] * B[j];
				}

				// calculate the reflection coefficient
				K = -(Inner / m_Diagonal[i-1]);

				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

               	// update predictor coefficients
				A[0] = ComplexFloat.Zero;
				Array.Copy(B, 0, A, 1, i);
				for (j = 0, l = i - 1; j < i; j++, l--)
				{
					A[j] += K * B[l];
				}

				// update prediction error
				m_Diagonal[i] = (ComplexFloat.One - K * K) * m_Diagonal[i - 1];

				// check if prediction error is zero (ie singular sub-matrix)
				if (m_Diagonal[i] == ComplexFloat.Zero)
				{
					m_IsSingular = true;
					return;
				}

			}

			// solution completed
			m_IsSingular = false;
			return;

		}

		#endregion Protected Members

        #region Public Methods
		
		/// <summary>
		/// Get a copy of the Toeplitz vector.
		/// </summary>
		/// <remarks>
		/// This member returns a copy of the left-most column of the Toeplitz matrix.
		/// </remarks>
		public ComplexFloatVector GetVector()
		{
			return new ComplexFloatVector(m_ToeplitzColumn);
		}

		/// <summary>
		/// Get a copy of the Toeplitz matrix.
		/// </summary>
		/// <remarks>
		/// Based on the internal vector this member constructs a square Toeplitz matrix.
		/// </remarks>
		public ComplexFloatMatrix GetMatrix()
		{
			int i, j;
				
			// allocate memory for the matrix
			ComplexFloatMatrix tm = new ComplexFloatMatrix(m_Order);

#if MANAGED
            // fill top row
			ComplexFloat[] top = tm.data[0];
			Array.Copy(m_ToeplitzColumn.data, 0, top, 0, m_Order);

			if (m_Order > 1)
			{
				// fill bottom row (reverse order)
				ComplexFloat[] bottom = tm.data[m_Order - 1];
					
				for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
				{
					bottom[i] = m_ToeplitzColumn[j];
				}
					
				// fill rows in-between
				for (i = 1, j = m_Order - 1 ; j > 1; i++)
				{
					Array.Copy(top, 0, tm.data[i], i, j--);
					Array.Copy(bottom, j, tm.data[i], 0, i);
				}
			}
#else     
			if (m_Order > 1) {
				ComplexFloat[] top = new ComplexFloat[m_Order];
				Array.Copy(m_ToeplitzColumn.data, 0, top, 0, m_Order);
				tm.SetRow(0, top);
                
				// fill bottom row (reverse order)
				ComplexFloat[] bottom = new ComplexFloat[m_Order];

				for (i = 0, j = m_Order - 1; i < m_Order; i++, j--) {
					bottom[i] = m_ToeplitzColumn[j];
				}

				// fill rows in-between
				for (i = 1, j = m_Order - 1 ; j > 0; i++){
					ComplexFloat[] temp = new ComplexFloat[m_Order];
					Array.Copy(top, 0, temp, i, j--);
					Array.Copy(bottom, j, temp, 0, i);
					tm.SetRow(i, temp);
				}
			} else {
				Array.Copy(m_ToeplitzColumn.data, 0, tm.data, 0, m_Order);
			}
#endif
			return tm; 				
		}
		
		///<summary>
		/// Get the lower triangle matrix of the UDL factorisation.</summary>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		public ComplexFloatMatrix GetLower()
		{
			// make sure the factorisation is completed
			Compute();
			
			// check if there was a result
			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}

			// copy jagged array into a ComplexFloatMatrix
			ComplexFloatMatrix Lower = new ComplexFloatMatrix(m_Order);
            for (int i = 0; i < m_Order; i++ )
            {
#if MANAGED
                m_LowerTriangle[i].CopyTo(Lower.data[i], 0);
#else
                for( int j = 0; j < i+1; j++ ){
                    Lower.data[j*m_Order+i] = m_LowerTriangle[i][j];
                }
#endif
            }
			
			return Lower;
		}
		
		///<summary>
		/// Get the upper triangle matrix of the UDL factorisation.</summary>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		public ComplexFloatMatrix GetUpper()
		{
			// make sure the factorisation is completed
			Compute();
			
			// check if there was a result
			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}
			
			// copy jagged array into a ComplexFloatMatrix and then calculate the transpose
			ComplexFloatMatrix Lower = GetLower();
			
			return Lower.GetTranspose();
		}

		///<summary>
		/// Get the diagonal matrix of the UDL factorisation.</summary>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		public ComplexFloatMatrix GetDiagonal()
		{
			// make sure the factorisation is completed
			Compute();
			
			// check if there was a result
			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}

			// copy diagonal vector into a ComplexFloatMatrix
			ComplexFloatMatrix Diagonal = new ComplexFloatMatrix(m_Order);
            for (int i = 0; i < m_Order; i++ )
            {
#if MANAGED
                Diagonal.data[i][i] = ComplexFloat.One / m_Diagonal[i];
#else
                Diagonal.data[i*m_Order+i] = ComplexFloat.One / m_Diagonal[i];
#endif
            }

			return Diagonal;
		}

		/// <summary>
		/// Get the determinant of the Toeplitz matrix.
		/// </summary>
		public ComplexFloat GetDeterminant()
		{
			ComplexFloat Determinant;

			// make sure the factorisation is completed
			Compute();
			
			// if any of the matrices is singular - just set the determinat to zero
			if (m_IsSingular == true)
			{
				Determinant = ComplexFloat.Zero;
			}
			else
			{
				Determinant = ComplexFloat.One;
				for (int i = 0; i < m_Order; i++)
				{
					Determinant *= 	m_Diagonal[i];
				}
			}

			return Determinant;
		}

        /// <overloads>
        /// Solve a symmetric square Toeplitz system.
        /// </overloads>
		/// <summary>
		/// Solve a symmetric square Toeplitz system with right-side vector.
		/// </summary>
		/// <param name="Y">The right-hand side of the system.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// Parameter <B>Y</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of Y is not equal to the number of rows in the Toeplitz matrix.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
		/// a symmetric square Toeplitz matrix, <B>X</B> is the unknown solution vector
		/// and <B>Y</B> is a known vector.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation
		/// using the Levinson algorithm, before calculating the solution vector.
		/// </para>
		/// </remarks>
        public ComplexFloatVector Solve(ComplexFloatVector Y)
        {

			ComplexFloatVector X;

            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y");
            }
            else if (m_Order != Y.Length)
            {
				throw new RankException("The length of Y is not equal to the number of rows in the Toeplitz matrix.");
            }

			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}

			int i, j, l;			// index/loop variables
			ComplexFloat Inner;			// inner product
			ComplexFloat G;				// scaling constant
			ComplexFloat[] A;				// reference to current order coefficients

			// allocate memory for solution
			X = new ComplexFloatVector(m_Order);

			// setup zero order solution
			X[0] = Y[0] / m_ToeplitzColumn[0];

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++)
			{
				// calculate inner product
           		Inner = Y[i];
           		for (j = 0, l = i; j < i; j++, l--)
           		{
           			Inner -= X[j] * m_ToeplitzColumn[l];
           		}
				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

				// update the solution vector
				G = Inner / m_Diagonal[i];
           		for (j = 0; j <= i; j++)
           		{
           			X[j] += G * A[j];
				}

			}

			return X;
        }

		/// <summary>
		/// Solve a symmetric square Toeplitz system with a right-side matrix.
		/// </summary>
		/// <param name="Y">The right-hand side of the system.</param>
		/// <returns>The solution matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// Parameter <B>Y</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The number of rows in Y is not equal to the number of rows in the Toeplitz matrix.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member solves the linear system <B>TX</B> = <B>Y</B>, where <B>T</B> is
		/// a symmetric square Toeplitz matrix, <B>X</B> is the unknown solution matrix
		/// and <B>Y</B> is a known matrix.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation
		/// using the Levinson algorithm, before calculating the solution vector.
		/// </para>
		/// </remarks>
        public ComplexFloatMatrix Solve(ComplexFloatMatrix Y)
        {

			ComplexFloatMatrix X;

            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y");
            }
            else if (m_Order != Y.ColumnLength)
            {
				throw new RankException("The numer of rows in Y is not equal to the number of rows in the Toeplitz matrix.");
            }

			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}

			int M = Y.RowLength;
			int i, j, l, m;			// index/loop variables
			ComplexFloat[] Inner;			// inner product
			ComplexFloat[] G;				// scaling constant
			ComplexFloat[] A;				// reference to current order coefficients
			ComplexFloat scalar;

			// allocate memory for solution
			X = new ComplexFloatMatrix(m_Order, M);
			Inner = new ComplexFloat[M];
			G = new ComplexFloat[M];

			// setup zero order solution
			scalar = ComplexFloat.One / m_ToeplitzColumn[0];
			for (m = 0; m < M; m++)
			{
#if MANAGED
                X.data[0][m] = scalar * Y.data[0][m];
#else
                X.data[m*m_Order] = scalar * Y.data[m*m_Order];
#endif
			}

			// solve systems of increasing order
			for (i = 1; i < m_Order; i++)
			{
				// calculate inner product
				for (m = 0; m < M; m++)
				{
#if MANAGED
                    Inner[m] = Y.data[i][m];
#else
                    Inner[m] = Y.data[m*m_Order+i];
#endif
				}

           		for (j = 0, l = i; j < i; j++, l--)
           		{
					scalar = m_ToeplitzColumn[l];
					for (m = 0; m < M; m++)
					{
#if MANAGED
                        Inner[m] -= scalar * X.data[j][m];
#else
                        Inner[m] -= scalar * X.data[m*m_Order+j];
#endif
					}
           		}

				// get the current predictor coefficients row
				A = m_LowerTriangle[i];

				// update the solution matrix
				scalar = ComplexFloat.One / m_Diagonal[i];
				for (m = 0; m < M; m++)
				{
					G[m] = scalar * Inner[m];
				}
           		for (j = 0; j <= i; j++)
           		{
					scalar = A[j];
					for (m = 0; m < M; m++)
					{
#if MANAGED
                        X.data[j][m] += scalar * G[m];
#else
                        X.data[m*m_Order+j] += scalar * G[m];
#endif
					}
				}
			}

			return X;
        }

		/// <summary>
		/// Get the inverse of the Toeplitz matrix.
		/// </summary>
		/// <returns>The inverted matrix.</returns>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices is singular.
		/// </exception>
		/// <remarks>
		/// This member inverts the matrix <B>T</B>,
		/// where <B>T</B> is symmetric square Toeplitz matrix.
		/// <para>
		/// The class implicitly decomposes the inverse Toeplitz matrix into a UDL factorisation
		/// using the Levinson algorithm, before using Trench's algorithm to complete
		/// the inversion.
		/// </para>
		/// </remarks>
        public ComplexFloatMatrix GetInverse()
        {
			Compute();

			if (m_IsSingular == true)
			{
				throw new SingularMatrixException("One of the principal sub-matrices is singular.");
			}

			ComplexFloatMatrix X;						// the solution matrix
			ComplexFloat[] A;							// prediction coefficients
			ComplexFloat Inner, gamma;
			int i, j, l, k, m, n;
			ComplexFloat[] XPrev;
			ComplexFloat AI;
			ComplexFloat AP;

			// allocate memory for the solution
			X = new ComplexFloatMatrix(m_Order);

			// get the current predictor coefficients row
			A = m_LowerTriangle[m_Order - 1];

			// calculate gamma
			Inner = ComplexFloat.Zero;
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
				Inner += m_ToeplitzColumn[i] * A[j];
			}
			gamma = ComplexFloat.One / Inner;

			// calculate first row of inverse
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
#if MANAGED
                X.data[0][i] = gamma * A[j];
#else
                X.data[i*m_Order] = gamma * A[j];
#endif
			}

            // calculate successive rows of upper wedge
			XPrev = X.GetRow(0).data;
			for (i = 0; i < m_Order / 2; i++, XPrev = X.GetRow(i+1).data)
			{
				AI = A[i];
				l = m_Order-2-i;
				AP = A[l];

				for (k = i + 1, m = i, n = 0; k < m_Order - i - 1; k++, l--, m++, n++ )
				{
#if MANAGED
                    X.data[i+1][k] = XPrev[m] + gamma * (AP * A[l] - AI * A[n]);
#else
                    X.data[k*m_Order+i+1] = XPrev[m] + gamma * (AP * A[l] - AI * A[n]);
#endif
				}
			}

			// this is symmetric matrix ...
			for (i = 0; i <= m_Order / 2; i++)
			{
				for (j = i + 1; j < m_Order - i; j++)
				{
#if MANAGED
                    X.data[j][i] = X.data[i][j];
#else
                    X.data[i*m_Order+j] = X.data[j*m_Order+i];
#endif
				}
			}

			// and a persymmetric matrix.
			for (i = 0, j = m_Order - 1; i < m_Order; i++, j--)
			{
				for (k = 0, l = m_Order - 1; k < j; k++, l--)
				{
#if MANAGED
                    X.data[l][j] = X.data[i][k];
#else
                    X.data[j*m_Order+l] = X.data[k*m_Order+i];
#endif
				}
			}

			return X;
        }

        #endregion Public Methods

        #region Public Static Methods

        /// <overloads>
        /// Solve a symmetric square Toeplitz system.
        /// </overloads>
        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side vector.
        /// </summary>
        /// <param name="T">The left-most column of the Toeplitz matrix.</param>
        /// <param name="Y">The right-side vector of the system.</param>
        /// <returns>The solution vector.</returns>
        /// <exception cref="ArgumentNullException">
        /// <B>T</B> and/or <B>Y</B> are null references
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <B>T</B> does not match the length of <B>Y</B>.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// A principal sub-matrix of <B>T</B> is singular.
        /// </exception>
        /// <remarks>
        /// This method solves the linear system <B>AX</B> = <B>Y</B>. Where
        /// <B>T</B> is a symmetric square Toeplitz matrix, <B>X</B> is an unknown
        /// vector and <B>Y</B> is a known vector.
        /// <para>
        /// The classic Levinson algorithm is used to solve the system. The algorithm
        /// assumes that all the principal sub-matrices of the Toeplitz matrix are
        /// non-singular. When a sub-matrix is near singular, accuracy will
        /// be degraded. This member requires approximately 2 times the order squared
        /// FLOPS to calculate a solution.
        /// </para>
        /// <para>
        /// See, <EM>G. H. Golub and C. F. Van Loan</EM>, Matrix Computations, Third
	    /// Edition, 1996, pages 196-197, The John Hopkins University Press.
        /// </para>
        /// </remarks>
        public static ComplexFloatVector Solve(ComplexFloatVector T, ComplexFloatVector Y)
        {

            ComplexFloatVector X;

            // check parameters
            if (T == null)
            {
                throw new System.ArgumentNullException("T");
            }
            else if (Y == null)
            {
                throw new System.ArgumentNullException("Y");
            }
            else if (T.Length != Y.Length)
            {
				throw new RankException("The length of T and Y are not equal.");
            }
            else
            {

                // allocate memory
                int N = T.Length;
                X = new ComplexFloatVector(N);                    // solution vector
                ComplexFloatVector a = new ComplexFloatVector(N - 1);   // prediction coefficients
                ComplexFloatVector Z = new ComplexFloatVector(N - 1);   // temporary storage vector
                ComplexFloat e;                                   // prediction error
                ComplexFloat g;                                   // reflection coefficient
                ComplexFloat inner;                               // inner product
                ComplexFloat k;
                int i, j, l;

                // setup zero order solution
                e = T[0];
                if (e == ComplexFloat.Zero)
                {
                    throw new SingularMatrixException("A principal sub-matrix is singular.");
                }
                X[0] = Y[0] / T[0];

                // calculate solution for successive orders
                for (i = 1; i < N; i++)
                {

                    // calculate first inner product
                    inner = T[i];
                    for (j = 0, l = i - 1; j < i - 1; j++, l--)
                    {
                        inner += a[j] * T[l];
                    }

                    // update predictor coefficients
                    g = -(inner / e);
                    for (j = 0, l = i - 2; j < i - 1; j++, l--)
                    {
                        Z[j] = a[j] + g * a[l];
                    }

                    // copy vector
                    for (j = 0; j < i - 1; j++)
                    {
                        a[j] = Z[j];
                    }

                    a[i - 1] = g;

                    e *= (ComplexFloat.One - g * g);
                    if (e == ComplexFloat.Zero)
                    {
                        throw new SingularMatrixException("A principal sub-matrix is singular.");
                    }

                    // calculate second inner product
                    inner = Y[i];
                    for (j = 0, l = i; j < i; j++, l--)
                    {
                        inner -= X[j] * T[l];
                    }

                    // update solution vector
                    k = inner / e;
                    for (j = 0, l = i - 1; j < i; j++, l--)
                    {
                        X[j] = X[j] + k * a[l];
                    }
                    X[j] = k;

                }

            }

            return X;

        }

        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side matrix.
        /// </summary>
        /// <param name="T">The left-most column of the Toeplitz matrix.</param>
        /// <param name="Y">The right-side matrix of the system.</param>
        /// <returns>The solution matrix.</returns>
        /// <exception cref="ArgumentNullException">
        /// <B>T</B> and/or <B>Y</B> are null references
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <B>T</B> does not match the number of rows in <B>Y</B>.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// A principal sub-matrix of <B>T</B> is singular.
        /// </exception>
        /// <remarks>
        /// This method solves the linear system <B>AX</B> = <B>Y</B>. Where
        /// <B>T</B> is a symmetric square Toeplitz matrix, <B>X</B> is an unknown
        /// matrix and <B>Y</B> is a known matrix.
        /// <para>
        /// The classic Levinson algorithm is used to solve the system. The algorithm
        /// assumes that all the principal sub-matrices of the Toeplitz matrix are
        /// non-singular. When a sub-matrix is near singular, accuracy will
        /// be degraded. This member requires approximately <B>N(M+1)N</B> FLOPS to calculate
        /// a solution, with <B>N</B> the order of the Toeplitz matrix and <B>M</B> the number
        /// of columns in the matrix <B>Y</B>.
        /// </para>
        /// <para>
        /// See, <EM>G. H. Golub and C. F. Van Loan</EM>, Matrix Computations, Third
	    /// Edition, 1996, pages 196-197, The John Hopkins University Press.
        /// </para>
        /// </remarks>
        /// <preliminary>
        /// This is an preliminary version. The speed can be improved, by paritioning
        /// the algorithm into two separate phases.
        /// </preliminary>
        public static ComplexFloatMatrix Solve(ComplexFloatVector T, ComplexFloatMatrix Y)
        {

            ComplexFloatMatrix X;

            // check parameters
            if (T == null)
            {
                throw new System.ArgumentNullException("T");
            }
            else if (Y == null)
            {
                throw new System.ArgumentNullException("Y");
            }
            else if (T.Length != Y.ColumnLength)
            {
				throw new RankException("The length of T and Y are not equal.");
            }
            else
            {

                // allocate memory
                int N = T.Length;
                int M = Y.RowLength;
                X = new ComplexFloatMatrix(N, M);                 // solution matrix
                ComplexFloatVector a = new ComplexFloatVector(N - 1);   // prediction coefficients
                ComplexFloatVector Z = new ComplexFloatVector(N);       // temporary storage vector
                ComplexFloat e;                                   // prediction error
                ComplexFloat p;                                   // reflection coefficient
                ComplexFloat inner;                               // inner product
                ComplexFloat k;
                int i, j, l, m;

                // setup zero order solution
                e = T[0];
                if (e == ComplexFloat.Zero)
                {
                    throw new SingularMatrixException("A principal sub-matrix is singular.");
                }
                for (m = 0; m < M; m++)
                {
                    X[0, m] = Y[0,m] / T[0];
                }

                // calculate solution for successive orders
                for (i = 1; i < N; i++)
                {

                    // calculate first inner product
                    inner = T[i];
                    for (j = 0, l = i - 1; j < i - 1; j++, l--)
                    {
                        inner += a[j] * T[l];
                    }

                    // update predictor coefficients
                    p = -(inner / e);
                    for (j = 0, l = i - 2; j < i - 1; j++, l--)
                    {
                        Z[j] = a[j] + p * a[l];
                    }

                    // copy vector
                    for (j = 0; j < i - 1; j++)
                    {
                        a[j] = Z[j];
                    }

                    a[i - 1] = p;
                    e *= (ComplexFloat.One - p * p);

                    if (e == ComplexFloat.Zero)
                    {
                        throw new SingularMatrixException("A principal sub-matrix is singular.");
                    }

                    // update the solution matrix
                    for (m = 0; m < M; m++)
                    {

                        // retrieve a copy of solution column
                        for (j = 0; j < i; j++)
                        {
                            Z[j] = X[j, m];
                        }

                        // calculate second inner product
                        inner = Y[i, m];
                        for (j = 0, l = i; j < i; j++, l--)
                        {
                            inner -= Z[j] * T[l];
                        }

                        // update solution vector
                        k = inner / e;
                        for (j = 0, l = i - 1; j < i; j++, l--)
                        {
                            Z[j] = Z[j] + k * a[l];
                        }
                        Z[j] = k;

                        // store solution column in matrix
                        for (j = 0; j <= i; j++)
                        {
                            X[j, m] = Z[j];
                        }

                    }

                }

            }

            return X;

        }

		/// <summary>
		/// Solve the Yule-Walker equations for a symmetric square Toeplitz system
		/// </summary>
		/// <param name="R">The left-most column of the Toeplitz matrix.</param>
		/// <returns>The solution vector.</returns>
		/// <exception cref="ArgumentNullException">
		/// R is a null reference.
		/// </exception>
		/// <exception cref="ArgumentOutOfRangeException">
		/// The length of R must be greater than 1.
		/// </exception>
        /// <exception cref="SingularMatrixException">
        /// A principal sub-matrix of <B>R</B> is singular.
        /// </exception>
		/// <remarks>
		/// This member is used to solve the Yule-Walker system <B>AX</B> = -<B>a</B>,
		/// where <B>A</B> is a symmetric square Toeplitz matrix, constructed
		/// from the elements <B>R[0]</B>, ..., <B>R[N-2]</B> and
		/// the vector <B>a</B> is constructed from the elements
		/// <B>R[1]</B>, ..., <B>R[N-1]</B>.
		/// <para>
		/// Durbin's algorithm is used to solve the equation. It requires
		/// approximately the order squared <EM>FLOPS</EM> to calculate the
		/// solution. The Yule-Walker system is used for linear prediction analysis
		/// and by the Trench algorithm for matrix inversion.
		/// </para>
		/// <para>
		/// See, <EM>G. H. Golub and C. F. Van Loan</EM>, Matrix Computations, Third
		/// Edition, 1996, pages 194-196, The John Hopkins University Press.
		/// </para>
		/// </remarks>
        public static ComplexFloatVector YuleWalker(ComplexFloatVector R)
        {

            ComplexFloatVector a;

            // check parameters
            if (R == null)
            {
                throw new System.ArgumentNullException("R");
            }
            else if (R.Length < 2)
            {
                throw new System.ArgumentOutOfRangeException("The length of R must be greater than 1.");
            }
            else
            {

                int N = R.Length - 1;
                a = new ComplexFloatVector(N);                    // prediction coefficients
                ComplexFloatVector Z = new ComplexFloatVector(N);		// temporary storage vector
                ComplexFloat e;				                    // predictor error
                ComplexFloat inner;                               // inner product
                ComplexFloat g;                                   // reflection coefficient
                int i, j, l;

                // setup first order solution
                e = R[0];
                if (e == ComplexFloat.Zero)
                {
                    throw new SingularMatrixException("One of the principal sub-matrices is singular.");
                }
                g = -R[1] / R[0];
                a[0] = g;

                // calculate solution for successive orders
                for (i = 1; i < N; i++)
                {

                    e *= (ComplexFloat.One - g * g);
                    if (e == ComplexFloat.Zero)
                    {
                        throw new SingularMatrixException("One of the principal sub-matrices is singular.");
                    }

                    // calculate inner product
                    inner = R[i + 1];
                    for (j = 0, l = i; j < i; j++, l--)
                    {
                        inner += a[j] * R[l];
                    }

                    // update prediction coefficients
                    g = -(inner / e);
                    for (j = 0, l = i - 1; j < i; j++, l--)
                    {
                        Z[j] = a[j] + g * a[l];
                    }

                    // copy vector
                    for (j = 0; j < i; j++)
                    {
                        a[j] = Z[j];
                    }

                    a[i] = g;

                }
            }

            return a;

        }

		/// <summary>
		/// Invert a symmetric square Toeplitz matrix.
		/// </summary>
		/// <param name="T">The left-most column of the symmetric Toeplitz matrix.</param>
		/// <returns>The inverted matrix.</returns>
		/// <exception cref="ArgumentNullException">
		/// <B>T</B> is a null reference.
		/// </exception>
		/// <exception cref="RankException">
		/// The length of <B>T</B> must be greater than 1.
		/// </exception>
		/// <exception cref="SingularMatrixException">
		/// One of the principal sub-matrices of T is singular.
		/// </exception>
		/// <remarks>
		/// This member inverts the matrix <B>T</B>,
		/// where <B>T</B> is symmetric square Toeplitz matrix.
		/// <para>
		/// Trench's algorithm is used to perform the inversion. It requires
		/// approximately 3.25 times the order squared <EM>FLOPS</EM> to calculate the
		/// solution. This algorithm is not suitable for systems with a matrix
		/// which has a singular or near-singular principal sub-matrix.
		/// </para>
		/// <para>
		/// See, <EM>G. H. Golub and C. F. Van Loan</EM>, Matrix Computations, Third
		/// Edition, 1996, pages 197-200, The John Hopkins University Press.
		/// </para>
		/// </remarks>
        public static ComplexFloatMatrix Inverse(ComplexFloatVector T)
        {

            ComplexFloatMatrix X;

            // check parameters
            if (T == null)
            {
                throw new System.ArgumentNullException("T");
            }
            else if (T.Length < 2)
            {
                throw new System.RankException("The length of T must be greater than 1.");
            }
            else
            {

				int N = T.Length;
				ComplexFloat f, g;
				int i, j, l, k, m, n;
				X = new ComplexFloatMatrix(N);

				// calculate the predictor coefficients
				ComplexFloatVector Y = ComplexFloatSymmetricLevinson.YuleWalker(T);

				// calculate gamma
				f = T[0];
				for (i = 1, j = 0; i < N; i++, j++)
				{
					f += T[i] * Y[j];
				}
				g = ComplexFloat.One / f;

				// calculate first row of inverse
				X[0, 0] = g;
				for (i = 1, j = 0; i < N; i++, j++)
				{
					X[0, i] = g * Y[j];
				}

				// calculate successive rows of upper wedge
				for (i = 0, j = 1, k = N - 2; i < N / 2; i++, j++, k--)
				{
					for (l = j, m = i, n = N-1-j; l < N - j; l++, m++, n--)
					{
						X[j, l] = X[i, m] + g * (Y[i] * Y[m] - Y[k] * Y[n]);
					}
				}

				// this is symmetric matrix ...
				for (i = 0; i <= N / 2; i++)
				{
					for (j = i + 1; j < N - i; j++)
					{
						X[j, i] = X[i, j];
					}
				}

				// and a persymmetric matrix.
				for (i = 0, j = N - 1; i < N; i++, j--)
				{
					for (k = 0, l = N - 1; k < j; k++, l--)
					{
						X[l, j] = X[i, k];
					}
				}

			}

			return X;

        }

        #endregion Public Static Methods

    }

}
